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PART  II:  THE  SYMMETRIC  POSITIVE  UEFINIIE  PROBLEM 
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requirements  as  well  as  Che  computational  effort  are  strongly  affected 


unknoms  are  ellaioaced.  A difficulty  in  using  iterative  methods  is  that  the  behavior  of  the 


iterative  methods  is  comparatively  recent 


behavior  than  methods  like  SOR.  A final  attractive  feature  is  that  the  model  problem.  It  arises  from  the  simplest  finite  difference 

2— norm  of  the  error  vector  decreases  monotonically  at  each  iteration  approximations  to  Poisson's  equation  in  a p— dimensional  cube  and  has 


elliptic  partial  differential  equations.  We  call  this  the  generalized  ^ computational  version  of  these  methods  and  discuss  tw  important 

model  problem”.  The  second  example  is  a special  case  of  the  generalized  special  cases  - the  Conjugate  Gradient  (OG)  and  the  Conjugate  Residual 


(CR)  methods.  In  Chapter  5,  we  motivate  and  apply  the  Idea  of  using  and  PCC  methods  with  some  competing  direct  and  Iterative  methods. 


residual  over  subspaces  of  lacreasia^  dimension.  It  can  be  viewed  as 


consistent  if  the  following  holds:  If  for  some  1,  x.  Is  a solution,  say 


Integer  s > 0,  ^ Is  Independent  of  1 for  all 


U(1  , 1,,...,  1 ) which  Is  an  approximation  to  ud.h,  l.h 


smooth  functions  with  o (x)  > m for  some  positive  constant 


+ h/2),  The  resulting  matrix  A is  block-trldlagonal 
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Figure  2.1:  Natural  ordering  for  points  on  a 3 x 3 mesh  on  the  unit  “®  obtain  the  symmetric  poaltlve-def Inlte  model  problem  In  two 

square  ( n-3,  N-9,  h«l/4  ) 

dimensions.  For  o sufficiently  positive,  A Is  Indefinite  and  we  obtain 


Is  ic(A)  « 0(n  ) [62].  For  the  Indefinite  model  problem,  the  eigenvalues 


where  >”  are  Che  eigenvalues  of  the  aatrix  corresponding  to  the  Iterative  methods  characteristically  require  an  Initial  guess 


n)  multiplications  are  required  to  compute  a matrlx-vector 


product.  If  A l8  positive  definite,  we  usually  use  a diagonal  scaling 


of  increasing  dimension.  These  methods  are  optimal  with  respect  to  the 


error  functional,  anong  all  linear  Iterative  methods.  Corresponding  to 


direction  vectors,  the  uni-dlrectlonal  minimizations 


using  the  Induction  hypothesis,  this 


and  hence 


and  Ap  Is  orthogonal  to  A p.»  j ^ i*  by  Theorem 


The  fundanental  minlBlzatlon  property  of  Theorem  3.3  allows 


(A))r, 


■in  ( £ t R (X  )v  , £ t.R  (X  ) X **  v ) ye  then  have  the  following  result  (see  (27]  for  the  case 


and  proceeding  as  In  the  proof  of  Theorem  3.5, 


we  have 


the  simplifications  to  the  general  algorithm  of  Chapter  3 due  to  the 


conjugate  residual  methods,  are  special  cases  of  this  class  of  methods* 


for  p > 0*  The  corresponding  number  of  multiplications  is 


Proof ; By  Theorea 


h-l/(n+l)  la  the  aesh-wldth  (341.  Thus  for  i < k.  Then  by  (3.10a)  and  (3.6d) 


A ^ Algorithm  4,1:  The  Variational  Method  for  Positive  Definite 


Stlefel  I39j  which  can  be  obtained  from 


Algorlttw  4,2: The  Conjugate  Gradient  Method  practical  point  of  view  as  Algorithm  4.2  since  It  requires  an  additional 


Rutlshauser  [27]  considered  a version  of  the  conjugate  graalent 


Each  Iteration  of  Algorithm  4,3  requires  6N  + 2 multiplications  If  the  measure  of  the  error  In  the  iteration  process  is  based  on 


(4,4)  II X - II  < II X - x^*^  II  . three-dimensional  model  problems  respectively  for  various  mesh-widths 


required  by  CG  to  reduce  the  2-aorm  of  the  error  by  a factor  of  1/n 


Table  4,4;  Work  required  by  the  CG  and  CR  methods  to  reduce  the  2-non 
of  the  error  by  10  for  model  problem  in  three  dimensions 


(a)  Two  dlneostona 


Graphs  Illustrating  0(n  log  n)  convergence  result  for  the 
CG  method  on  the  model  probleo 


Chapter  A could  be  used  to  find  aa  approxiinatloQ  to  x' • from  which  an  require,  assuming  no  roundoff  error,  at  most  m iterations  to  converge  to 


If  R has  only  m distinct  eigenvalues,  then  K*  also  has  some  vector  u,  then  the  overall  %rork  per  iteration  may  go  up 


and  work  raqulreaeocs. 


following  change  of  variables: 


process;  By  a simple  analysis. 


or,  equivalently,  “e  have 


Hence  (3.ila)  holds  for  1 * k'fl 


or  recursively 


the  matr iX'vector  product 


In  Section  2,  we  discuss  diagonal  scaling,  which  is  equivalent  to 
choosing  K to  be  a diagonal  matrix.  In  Section  3,  we  give  examples  of 
problems  for  which  the  splitting  can  be  chosen  such  that  m”^A  has  only 


be  derived  for  Che  ocher  precondicloned  mechods  Incroduced 


In  Section  4,  we  develop  a particular  approach  of  "approximately"  In  Section  7,  we  present  the  results  of  numerical  experiments  which 

factoring  A for  the  generalized  model  problem.  Included  as  special  compare  the  performance  of  the  PCG  method  for  the  model  and  anisotropic 

cases  are  the  approximate  factor! Tations  of  Meljerink  and  van  der  Vorst  nroblems  for  the  different  nrecnnHitioninaa. 


compared  to  the  llne-SSOR  preconditioning  in  three  dimensions 


scaled  natrix  are  unity,  ve  save  N multiplications  each  time 


unit  upper  triangular,  and  L and  U have  non'-’zeroes  only  In  positions  in  Define 


diagonals  of  A, 


where  ve  assume  here  and  in  Che  following  sections  Chat  subscripted 


use  different  approximations* 


compute  an  approximate  factorization  LL  of  A by  the  Cholesky  algorithm 


factorization  procedure  proposed  by  Meijerink  and  structure  as  the  lo%^er  triangular  part 


0(d  ) multipllcatioas  In  three  dimensions,  we  have  the  following  result* 


asymptotically  the  same  as  the  condition  number  of  A.  where  a Is  a parameter.  With  this  definition  of  L and  U 


seI£*adJolnc  Neumann  problem  (with  slightly  different  conditions 


both  cases.  Hence  a result  similar  to  Theorem  6.3  can  be  stated  for 


5 


using  Stone's  preconditioning.  Numerical  results  showing  that  this 


stored,  and  each  solution  of  Mu  ■ v requires  backsolves  corresponding  to 


typical  situation,  H and  V would  be  trldlagonal,  or  could  be 


6.6:  The  BIock'SSOR  Precondiclonia 


so  that  M is  symmetric  and  positive  definite 


block-SSOR  method  is  a monotone  decreasing  function  of  P(S  ')»  where  Ue  now  consider  a few  particular  choices  of  the  partition  »•  Denote  by 


1+(1-2o(b'  ^)+4y) 


(b)  in  three  dimensions 


before  applying  Che  PCG  method  with  the  block^SSOR  preconditioning  gives 


In  general,  we  can  always  choose  the  scaling  matrix  as  in  (6.4)  so 


(6.23a)  A . W.  - V.  + w I a^.  clearly  requires  twice  the  number  of  multiplications  as  there  are 


For  Che  model  problem  we  do  doC  scale  the  system  as  we  want  to  take  Let 


diagonal.  There  are  some  overhead  operations,  but  these 


problem  where  Che  A are  all  different,  the  decomposition  has  to  be 


(371  that  makes  systems  with  them 


In  three  dlaenslons  plane  preconditioning  gives  subsystems  with  diagonals  are  equal  to 


Theorem  6»  7:  For  the  PCG  method  using  block-SSOR  preconditioning,  the 


be  greatly  reduced  by  going  from  point  to  line  and  line  to  plane 


the  two-'  and  three-dimensional  generalized  model  problems  respectively.  preconditioning  • the  total  work  required  is  greater.  However,  the 


three-dimensional  model  problems.  Tables  6.3a  and  6. 3b  give 


number  la  small  so  chat  Che  Iteration  process  shows  Improved 
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t The  plots  corresponding  to  the  two  preconditionings  were  alno 
Identical.  Therefore,  only  one  has  been  shown. 


112  - _ 113 


Figure  6« lb;  Error  reduction  for  the  CC  and  PCG  nethods  on  model  problem 
in  three  dimensions  (h  « 1/16) 


(a)  Dupont*  Kendall,  and  Rachford  Precondltlon:ing 


o.eoo  Q.eso  0.900  0.950  1.000  1.050 

Table  6.4;  Comparison  of  point-  and  llne-SSOR  preconditionings  for  the 
generalized  model  problem  (t»  * ^i)* 


(c)  Polnc-SSOR  PrecoDdltloalag 


le  6»6;  Comparison  of  ICCG(O)  and  1CCG(3)  for  model  problem 
in  two  dimensions 


Figure  6, 3;  Graphs  iliuscrating  convergence  behavior  of  the  PCG 
method  on  model  problem  in  two  dimensions 
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particular  preconditioning  of  the  two-cyclic  system.  In  Section 


that  it  is  equivalent  to  a method  due  to  Reid  [54]  in  the  sense  the  CG  iteration  requires  matrix-vector  products,  this  would  require 


wish  to  choose  a preconditioning  matrix  M =QQ  such  that 


t In  the  next  section  we  will  show  that  exactly  half  the  number  of 
iterations  are  required  if  CG  is  applied  to  (7.3)  rather  than  to  (7.2). 


storage  Is  requirec 


In  Algorithm  7.2,  Steps  2 and  3 actually  constitute  a double 


depend  on  the  for  their  computation.  Ue  now  show  that  there  is 


that  the  r'  can  be  generated  without  explicitly  forming 


Making  these  substitutions  in  Algorithm  7.2,  we  get  the  following.  Comparing  the  Algorithms  7.1  and  7.3,  it  is  clear  chat 


0 


Step  A;  Terminate  the  Iteration  process  and  compute  x, 


t The  work  requirements  for  the  RFCG  method  derived  in  Section  7.2  have 
been  halved,  since  the  RPCG  method  requires  exactly  half  the  number  of 
iterations  required  by  the  other  three  algorithms. 


Assume  witliout  loss  of  generality  Chat  S is  n,  x a,  and 


the  following  form  in  which  we  are  %rorkiag  directly  with  the  variables 


Uork/Xteratlon  Storage 


t The  work  requirements  for  reduced  CG  with  two-cycllc  preconditioning 
have  been  halved,  since  It  requires  exactly  half  the  number  of 
Iterations  required  by  the  other  two  algorithms. 


of  trldlagonal  matrices  In  ADI  and  the  Dupont,  Kendall,  and  Rachford 


The  number  of  Iterations  required  for  ADI  and  SOR  were  taken  from 
Birkiioff,  Varga,  and  Young  (5]  and  those  for  SSOR-SI  from  Ehrlich  [21]. 
For  SOR  and  SSOR-SI,  the  numbers  correspond  to  the  value  of  w that 
required  the  fewest  number  of  iterations*  For  PCG,  no  attempt  was  madi 
to  choose  a optimally;  a » 0 was  used  throughout* 


termlaated  Che  Iceracion  when  and  61  respectively.  For  this  choice  of  u,  PCG  Is  more  efficient  than 


should  oaerge  as  the  superior  method 


352K,  more  than  double  that  required  by  PCG. 


Che  chree-diaenslonal  model  problem  on  an  n x n x n mesh,  the 


•(«lton«ry  value  at  y iff 


If  A is  not  positive  definite*  then*  it  is  possible 


SYHMBK  algorithm  that  uses  this  factorization 
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Multiplying  out  the  left-hand  side  of  (9.20)  we  get 


\ 1 ®1  chosen  to  be  I x 1,  (9.22),  (9.24)  and  (9.25)  may 

0 be  used  to  compute  x from  x. 


Umlng  (9e22a)and  (9e24a),  we  get  from  (9s 19)  that 


(roB  now  on  we  suppress  Che  superscripts  on  the 


can  thus  say  that  each  Iteration  of  SYMMBK  requires  essentially 


For  Che  second  equality,  we  note  that  from  (9.24)  and  (9*28)  ve  can 


decomposition  to  solve  the  equations  (9.13).  SYMMLQ  was  proposed  by 

Since  D has  diagonal  blocks  of  order  1 or  2 only,  n , . 

1 Paige  and  Saunders  [52).  In  this  section,  we  briefly  outline  the  method 


Aod  show  that  Ic  is  less  efficleac  than  SYMMBK. 


of  roundoff  error,  when  x obtained  by  SYMMBK  it  can  be  obtained  from 


2M  fewer  nulcipllcacions  per  Iteratioa  as  conpared  to  SYMMLQ,  so  that  la 


CG  aethod  applied  to  <9.X)  breaks  down  if,  at  any  stage,  the 


* When  nonsingular,  an  iteration  of  Algorithm  9.3  requires 

^*^i+2*  ^i+1^ 

'i+2  " *’i+2 ^ p 6N  + 2 multiplications  plus  the  matrix  vector  product  Ap^,  When  p^^ 


properties  in  Theorem  9.3  will  not  hold  for  subsequent  iterations  and  it 


AO-A053  313 


UNCLASSIFIED 


YALE  UNIV  NEW  HAVEN  CONN  DEPT  OF  COMPUTER  SCIENCE  F/6  12/1 

CONJUGATE  GRADIENT  METHODS  FOR  PARTIAL  DIFFERENTIAL  EQUATIONS. (U) 

JAN  7S  R CHANDRA  N00014-76-C>0277 


RR-129 


NL 


2of2 

48m3I3 


solution.  The  algorltha  nay  even  have  serious  stability  problems 


either  nonsingular  or  the  second  nenber  of  a hyperbolic  pair 


^|^)/(Ap^,  pj^)  if  is  noDstaguiar 


showing  the  relation  between  the  direction  vectors  In  the 


9! 

a 


9 


1 


2 pivot  Is  chosen.  The  example  also  Illustrates  our  previous  result.  It  shows  Chat 


problaas.  In  Figure*  9.U  and  9.1b  ue  plot  the  error  versus  ounber  of  9.»b  that  for  a fixed  problem  (l.e.  fixed  o)  the  graph  of  the  number  of 


the  difference  approximations  are  only  second-order 


dinensions 
■ 1/32) 


Algorithm  4.2,  the  CG  method,  to  solve  the  normal  equations 


required  for  the  four  vectors  x,  r,  p and  Ap  • A^r  can  share  storage 


Applying  the  convergence  result  of  Theorem  4.1  to  (10.2)  with  We  now  apply  SYMMBK  (y  > 0)  to  the  linear  system  (10.3) 


the  first  m components  of  z,  and  by  z'  ' the  vector  consisting  of  the 


Lemma  10. 1:  If  i is  even,  o , i and  - 0.  If  l is  odd,  - 0 and 


k+1  “ ^'k+i 


(2).T 


We  now  prove  the  equivalence  of  the  two  aethoda 


k'i’K  We  do  this  by  considering  the  orthogonality  properties  of  the 


(10.8)  - f - A*^+1  - - AU.  where 


efficient  for  the  linear  least  squares  problem 


routine  application,  the  method  suffers  from  computational  difficulties 
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for  p > 0.  The  number  of  multlplicetlons 


Note  that  this  result  requires  that  the  two  Intervals 


obtain  bounds  on  the  rate  of  convergence. 


and  since  the  second  tens  on  the  right-hand  side  vanishes  Assune  (b)  and  (c)  hold  for  1 < k and 


by  the  Induction  hypothesis  we  have  that 


seep  1,  even  chough  Che  shorcer  one  Is  sufflclenc 


(A  Pj,p^) . Compuce 


(Ap^,A  P^)/(p|^»A  Pj^)*  Algorithm  IK  2:  The  Modified  Coalugate  Gradient  (HCR)  Method 


Step  4b;  Compute 


Compute  . (1^, Ap^) /(Ap^ .q^) 


would  like  to  choose  a preconditioning,  or  equivalently  the 


three-dimensional  indefinite  model  problems.  Values  of  a and  h were 


Tables  IU2  and  IK  3 give  Che  total  work  required  to  reduce  the  error  to 
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Figure  11.1;  Error  reduction  by  MCR  for  indefinite  model  probla 


three  dlnenslons  (a  • 100.  h • 1/16) 
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t For  the  model  problem  in  three  dimensions*  an  analogous  result  is 
obtained • 


Proof:  We  first  prove  (a)  by  Inductloo  on  1.  Since  Vj  a 4,  (a)  holds  We  now  prove  the  following  general  result  for  strictly  diagonally 


Then,  using  (a),  (c)  follows.  Then,  (see  Varga  (60J,  Theorem  3.7) 


Substituting  In  (A. 7),  the  result  follows 


the  condition  nimbcr  of  the  ICCG(O)  iteration  matrix  is  bounded  above 
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(36]  A.  van  der  Sluts*  Condition  numbers  and  equilibration  of 
matrices*  Numerische  Mathematlk  14:14-23,  1969. 


